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1.  Introduction 


In  the  prequel  (7)  to  the  present  paper  the  Fried  approach  to  short-exposure  imaging  through 
optical  turbulence  (2)  was  reanalyzed.  The  subject  of  that  study  was  the  short-exposure 
modulation  transfer  function  (MTF)  that  describes  the  mean  behavior  of  a  system  acquiring 
short-exposure  images  through  turbulence.  This  MTF  is  commonly  described  as  the  product  of  a 
system  MTF  and  an  atmospheric  MTF.  In  the  short-exposure  case,  the  net  MTF  is  expressed 
using  the  equation. 


Ms(u )  =  —  J  du  W\(u )  W\{u  —  u>)  x  exp  j  — 


\Di{D  |u;|)  + 
2 


where  u  and  uj  are  normalized  two-dimensional  vectors,  u  =  x/D  is  a  position  vector  in  a 
normalized  aperture,  while  u>  is  an  angular  frequency  vector  of  the  plane  perpendicular  to  the 
main  axis  of  propagation,  lu  =  |cj|  =  expresses  angular  frequency  relative  to  f20  =  D/X 

the  maximum  angular  frequency  (cyc/mrad)  passed  through  an  aperture  of  diameter  I)  by 
radiation  of  wavelength  A  (see,  e.g.,  Goodman  (3)).  W\  is  a  normalized  aperture  function, 


WD{x) 


1,  |x|  <  D/ 2; 

0,  \x\  >  D/2  ; 


(2) 


and  Vc(Duj)  and  D0(DD)  are  amplitude  and  phase  structure  functions. 


This  short-exposure  MTF  is  usually  written  as  the  product, 


Ms( uj)  =  M0(lu)  Msa(uj), 


(3) 


where,  M0( u)  is  the  system  MTF,  given  by, 


M0 (i 'jj )  =  —  J  du  iTj  (it)  ITj  (u  —  u))  = 


cos  1(cc)  —  CC\/1  —  <^2 

7/2 


W2(u ) 


The  key  contribution  in  the  former  paper  involved  the  retention  of  a  tilt-phase  correlation  term 
that  produced  turbulence  and  diffraction  related  effects  on  the  short-exposure  correction  to  the 
atmospheric  MTF.  The  present  paper  improves  upon  this  expression  by  exploring  further  into  the 
space  of  a  fourth-order  integral  introduced  in  the  previous  analysis.  Further  investigation  into  the 
behavior  of  this  integral  has  revealed  a  new  dependence. 
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To  describe  this  new  dependency,  we  first  consider  Fried’s  original  formulation  of  the 
short-exposure  atmospheric  MTF  (SAMTF).  According  to  the  standard  approach,  the  SAMTF  is 
itself  formed  as  a  correction  to  the  long-exposure  atmospheric  MTF  (LAMTF),  which  is  written, 

Mla(u)  =  exp{  — (2.1  A"  cu)5/3},  (5) 

wherein  the  factor  X  =  D  jrQ  represents  the  ratio  of  the  system  entrance  pupil  diameter  to  the 
Fried  coherence  length  (2),  which  for  homogeneous  Kolmogorov  turbulence  may  be  expressed  as, 

r0  =  2.1  x  1.437  (k2  LC^~V\  (6) 

with  k  =  2tt/\  the  optical  wavenumber,  L  the  path  length,  and  C'2  [m-2/3]  the  refractive  index 
structure  parameter.  From  the  LAMTF,  the  SAMTF  has  been  expressed  as  a  perturbation, 

Msa(oo)  =  exp{  — (2.1  X)5/3  [a;5/3  -  V(Q,  X)  a;2]},  (7) 

where  V (0.  X)  was  modeled  as  a  constant  in  Fried’s  original  analysis  (denoted  as  a),  taking 
values  of  either  1/2  or  1.  In  the  re-analysis  of  Tofsted  (7),  V ( Q ,  X)  was  shown  to  vary  with  both 
X  and  a  diffraction-related  parameter  Q  =  D/ P,  where  P  =  (A  L  ) 1/2  is  a  Fresnel  scaling 
parameter.  Together,  the  X  and  Q  parameters  account  for  turbulence  and  diffraction  influences 
on  V. 


The  meaning  of  this  short-exposure  correction  derives  from  Fried’s  observation  that 
long-exposure  imaging  averages  over  all  phase  aberrations  in  the  system  aperture,  but  when 
performing  short-exposure  imaging  the  tilt  of  the  incident  wave  may  shift  the  centroid  of  a  given 
object  point,  but  would  not  affect  its  image  clarity.  The  correction  factor,  V (Q.  X)  u2,  removes 
this  tilt  effect,  which  accounts  for  both  a  tilt- variance  effect  and  a  tilt-phase  correlation  term, 

V(Q,X)  =  N(Q,X )  -  1.0433  G(Q).  (8) 


where  N(Q,  X)  represented  the  tilt-phase  term,  and  G(Q)  represented  the  diffraction-influenced 
tilt- variance  effect. 


V(Q,  X)  was  approximated  in  the  prequel  using  the  analytic  form, 


V(Q,X)  «  A  + 


B 

10 


exp 


(x  +  l)3 
3.5 


(9) 


2 


where  x  =  log10(A"),  and  functions  A  and  B  were  modeled  as  functions  of  Q,  using  a 
log-variable  version,  q  =  log2(Q): 


A 


0.840  +  0.280  E(qa),  qa  =  0.51  (q  +  1.50);  q  <  -1.50 
0.840  +  0.116  E(gb),  qb  =  1.35  (q  +  1.50);  q  >  -1.50 


(10) 


B  =  0.805  +  0.265  E(gc),  qc  =  1.45  (q  m  0.15),  (11) 

where  E (z)  =  [exp(z)  —  l]/[exp(z)  +  1]  is  a  sigmoidal  function.  The  behavior  of  V(Q,  X )  is 
plotted  in  figure  1  which  is  a  reproduction  of  figure  4  of  Tofsted  (i). 
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Figure  1.  Intercomparison  of  V (Q,  X )  results  for  varying  Q  as  a 
function  of  X  from  the  prequel. 


Subsequent  to  the  publication  of  Tofsted  ( 1 ),  questions  arose  regarding  the  validity  of  the 
V(Q,X)  function  and  the  implication  that  super-resolution  results  were  obtained  (4).  Review 
and  reinvestigation  have  confirmed  elements  of  Charnotskii’s  critique  and  have  led  to  the  current 
paper. 

In  the  reinvestigation,  it  was  first  noted  that  equation  1  could  be  formed  through  a  combination  of 
equations  28,  29,  and  56  of  Tofsted  (/),  but  the  structure  functions  within  the  exponential  are 
positive  definite,  representing  mean  squared  quantities.  However,  if  V ( Q ,  A")  exceeds  unity  MSa 
will  exceed  unity  at  high  c u,  as  it  does  when  X  <  2  <  Q  (see  figure  1). 
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The  calculation  of  V (Q,  X)  primarily  relied  on  evaluation  of  N(Q,  X )  based  on  the  model, 


Mt(Q,  X,  u)  «  exp  [(2.1  X)5/3  N(Q,  X)  iJ 


(12) 


where  MT{Q,  X ,  to)  was  given  by  the  integral  relation  (equation  49  of  Tofsted  (/)), 


Mt(Q,  X,  u)  =  -  f  du  — ^  exp{  +  F(jD’  ^)j, 

7T  J  M0(u>)  I  I 


(13) 


32  r 

F(D,  it,  w)  =  —  /  duaWi(ua)  [wa  ■  cj] 

7 r  J 

x  (+£>4-°  (u-ua-  cj)]  -  £>4-°  (it  -  wa)]) .  (14) 

Due  to  circular  symmetry,  Mt  only  depends  on  the  magnitude  u  =  |iu|.  That  is,  one  may  always 
orient  lo  along  the  x-axis  without  loss  of  generality. 

In  the  prequel,  based  on  the  model  given  by  equation  12,  and  based  on  figure  3  of  Tofsted  ( 1 )  that 
appeared  to  show  cu2  dependence  of  the  exponent,  it  was  believed  to  be  unnecessary  to  evaluate 
equation  13  at  every  value  of  u  in  order  to  evaluate  a  single  point  V (0.  A").  Rather,  a  single 
evaluation  of  MT(Q,  X ,  uS)  could  be  performed  for  each  0  and  X  combination,  at  a 
“characteristic”  value  of  c o,  call  it  ujk.  N (Q,  X)  could  then  be  evaluated  for  this  computed  MT 
using  the  equation  12  model  and  the  value  of  u>k-  The  value  of  uk  was  chosen  so  that  it  would 
have  the  greatest  impact  on  the  overall  curve,  and  thus  the  point  Ms(u>k )  =  1/2  was  selected  to 
model  the  behavior  in  the  immediate  vicinity  of  the  knee  of  the  complete  MTF  curve. 

This  choice,  it  turns  out,  had  one  useful  feature,  and  one  detrimental  feature.  The  useful  feature 
is  that  while  (as  I  will  show)  we  must  abandon  the  previous  model  in  general,  as  an  approximation 
it  is  actually  rather  good.  When  compared  with  the  newer  model  described  further  along  in  the 
current  paper,  it  exhibits  only  a  third  of  the  error  of  the  original  Fried  equation.  The  detrimental 
feature  is  that  as  the  X  parameter  increases,  the  knee  of  the  Ms  curve  falls  at  ever  smaller  values 
of  u.  Thus  the  evaluation  point  is  at  smaller  ui,  and  if  (as  we  shall  find)  the  N  function  is  a 
varying  function  of  to,  then  the  correlation  between  uoK  and  X  may  produce  unusual  effects. 

A  further  reason  for  restricting  the  computations  to  an  evaluation  at  cuk  was  a  numerical  precision 
issue  with  evaluating  equation  14  and  its  effect  upon  introducing  this  result  in  the  exponential 
argument  in  equation  13.  The  phase  structure  function  is  modeled  as, 

V^Du)  =2(2.1X)5/3M5/3a(Q«),  (15) 
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where  a(Q  u )  is  a  diffraction  related  function  that  varies  smoothly  between  1/2  and  1  over  an 
8-decade  range  of  its  argument  (see  Tofsted  (7)).  However,  the  X 5/3  factor  can  produce 
exponential  arguments  exceeding  105  in  equation  13.  In  the  prequel  the  integration  process  was 
handled  on-the-fly,  which  we  found  to  be  computationally  intensive.  To  resolve  the  question  of 
how  N  varies  as  N(Q,  X,  u),  this  integration  routine  was  completely  reworked,  as  described  in 
section  2.  Then,  in  section  3  the  updated  V ( Q ,  X,  to)  function  is  considered,  and  its  behavior  is 
modeled  in  an  approximate  analytical  form,  as  previously,  but  now  as  a  function  of  three 
variables.  The  resulting  behavior  of  the  new  functional  form  is  considered  in  section  4. 


2.  Revised  Integration  Technique 


To  explore  the  variations  of  V (Q,  X,  u)  across  the  effective  spectrum  of  Q,  X,  and  u  we  consider 
Q  values  extending  from  1/1 6th  through  16  by  factors  of  2  (9  points),  X  from  0.1  to  103  in  steps 
of  1/1 0th  of  a  decade  (41  points),  and  u  from  0.01  through  0.99  plus  a  few  extra  points  at  small  cu 
(103  points).  In  sum,  the  required  number  of  calculated  points  increased  by  a  factor  of  over  100. 
Hence,  a  new  integration  technique  was  needed  that  was  both  faster  and  more  accurate. 

The  main  problem,  of  course,  was  the  factor  (2.1  A")5//3  =  S 5//3  in  D,,,  inside  the  F  double  integral 
of  equation  14.  To  deal  with  this  situation,  let  us  first  rewrite  equation  15  as, 

TV  Du)  =  2(2.1  X/Qf/z<p{Qu),  (16) 

where, 

<p(x)  =  x5/3  a(x).  (17) 

Using  this  functional  form  for  the  phase  structure  function,  equation  14  can  be  rewritten, 

32  r 

F(D,  m,  ut)  =  —  dua  Wi(ua)  [ua-oj] 

71  J 

x  i+V^D  (u-ua  —  u>)]  -  V,p[D  (u  -  Ma)]) 

=  ^(2.1  X/Qf/3J  dUaW^Ua)  [u0-«] 

X  (+<fi[Q  (u  -  Ua  -  u?)]  -  (p[Q  (u  -  Ma)]) 

=  (2.1  X)5/3  F(Q,  u,  cjx).  (18) 

This  first  step  allows  us  to  factor  out  (2.1  X)5/3  from  the  inner  double  integral.  The  function  F 
can  then  be  tabulated  for  |m|  <  1/2,  0<tc<l  (where  radial  symmetry  always  permits  us  to 
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orient  u>  along  the  +x  axis),  and  Q,  thereby  reducing  the  problem  of  evaluating  equation  13  from 
a  four-dimensional  integral  to  a  two-dimensional  integral  plus  an  interpolation  on  a 
multi-dimensional  look-up  table  of  F(...)  results.  This  greatly  increases  the  calculation  speed. 


The  remaining  major  problem  due  to  the  X5/3  factor  is  solved  by  the  addition  of  an  iterative  step 
during  the  computation  of  the  outer  double  integral.  The  problem  is  that  the  Mr  integral  must  be 
computed  with  (2.1  X)5^3  F(...)  inside  the  exponential  factor.  And  while  F  evaluates  to  a 
quantity  less  than  10,  the  A"5/3  factor  can  be  up  to  105.  And  yet,  once  Mr  is  evaluated,  the 
function  N  is  evaluated  by  taking  the  log  of  M> r  and  dividing  out  the  same  A"5/3  factor: 


N(Q,  X ,  u) 


In  [Mt(Q,  X,  u;)] 
(2.1  X)5/3  u2 


(19) 


Thus  the  (2.1  X)5/3  factor  appears  to  cancel,  but  cannot  simply  be  removed  from  the  outer  double 
integral  prior  to  actually  evaluating  Mr-  Nevertheless,  one  can  write  the  Mr  integral  in  the  form, 


. .  4  exp(+A)  r  ,  W\ (u  —  u>)  W\ (u)  ,  , 

Mt{Q,  X,  u)  =  -  /  du - — — - -  exp  +  F(D,  u,  u>)  -  A 

7T  J  Mq(U) 

=  Mjj(Q,  X,  u)  exp(+A), 


(20) 


such  that, 


N(Q,  X,  u) 


In  [Mu{Q,  X,  u)]  +  A 
(2.1  X)5/3  cu2 


(21) 


where  a  constant  such  as  exp(+A)  may  always  be  both  multiplied  and  divided  from  a  quantity. 
This  is  necessary  because  the  resolution  limit  of  double  precision  does  not  extend  to  exp(105),  yet 
a  priori  we  do  not  know  how  much  to  initially  subtract  inside  the  exponent.  And  when  we 
actually  try  to  evaluate  the  result  numerically,  due  to  the  double  precision  limitation  we  only 
know  that  the  result  has  gone  out  of  its  proper  bounds,  but  not  by  how  much.  We  must  therefore 
follow  an  iterative  procedure  to  find  an  appropriate  offset.  Thus  we  initially  set  A  =  0  and 
attempt  the  numerical  computation.  Each  time  the  algorithm  produces  an  out-of-bounds  result,  A 
is  incremented  by  +40,  and  the  numerical  integration  routine  is  run  again.  The  integration  is 
repeated,  adjusting  A,  until  such  point  that  A  is  adjusted  so  that  the  numerical  integration  result 
produces  a  meaningful  answer  in  the  double  precision  logic  of  the  routine. 


Given  both  the  tabulated  and  interpolated  inner  integral  approach,  along  with  the  outer  integral 
adjusted  by  the  offset  A,  a  complete  set  of  V (Q,  X,  uS)  values  could  be  computed.  Results  were 
computed  over  the  ranges  of  Q,  X,  and  uj  variations  indicated  at  the  beginning  of  this  section. 
This  range  of  the  parameters  covers  the  majority  of  the  dynamic  range  of  variations  for  most 
propagation  conditions  and  optical  systems  used  in  terrestrial  imaging  scenarios. 
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Figure  2  illustrates  the  two  extremes  of  the  data  set  by  plotting  two  sets  of  results  for  different  x 
values  as  functions  of  u  for  Q  =  1/16  and  16.  (Only  half  the  curves  are  labeled.)  From  these 
curves,  the  values  of  V(Q,  X ,  0)  appear  to  be  a  function  of  0  alone.  And  apparently 
V (Q,  X,  1)  ~  V (().  X,  0)  to  within  the  level  of  accuracy  of  the  calculations.  Also,  families  of 
curves  developed  for  different  Q  values  exhibit  similar  features,  indicating  that  further  analysis 
(in  the  following  section)  may  be  capable  of  modeling  this  behavior  using  fewer  than  the  38,007 
elements  in  the  calculated  sample  database.  The  result  of  said  analysis  is  a  fitting  function  using 
a  series  of  equations  that  produces  estimates  ofV{Q,  X ,  to)  that  are  virtually  identical  to  the 
computed  results. 


These  results  can  be  compared  with  the  behavior  of  the  previous  analysis  results  seen  in  figure  1 . 
A  critical  difference  is  that  unlike  the  previous  analysis,  V (Q,  X,  u)  does  not  degrade  with 
increasing  X  as  did  V (Q,  A").  The  reason  the  previous  analysis  showed  decreasing  V ((),  X) 
with  increasing  X  was  because  the  evaluation  point,  ujk,  was  a  decreasing  function  of  A". 


3.  Analysis  of  V (Q,  X ,  ui)  Parameterization 


The  analysis  process  whereby  the  V (Q,  X,  u)  dependence  was  determined  involved  a  lengthy 
sequence  of  analytical  results  to  unpack  the  general  behavior  of  the  function  into  a  sequence  of 
similar  behaviors.  The  first  stage  of  this  reduction  involved  the  division  of  the  function  into  two 
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segments:  the  base  segment,  designated  V0,  apparently  a  function  of  only  0.  and  Va,  a  variable 
portion  dependent  on  Q,  X,  and  u: 

V(Q,  X,  u)  =  V0(Q )  +  Va(Q,  X ,  u).  (22) 


To  model  several  of  the  resulting  functions  involved  in  characterizing  these  two  constituent 
functions  let  us  introduce  a  sigmoidal  function, 

[exp(g)  -  1]  =  [exp(+o;/2)  -  exp(-s/2)]  _  , 

[exp  (a;)  +  1]  [exp(+x/2)  +  exp(— x/2)\ 

We  then  employ  a  scaled  and  shifted  form  of  this  function: 

S(x,A,B,C,D)  =  A  +  BE[C(x-  D)\.  (24) 


We  shall  also  frequently  need  a  spliced  version  of  this  function  in  which  the  exponential  decays 
on  either  side  of  the  x  =  D  center  appear  to  exhibit  different  forms,  designated, 


Ss{x,A,BuB2,CuC2,D) 


A  +  Bi  £[Ci  (x  —  D)\,  x  <  D; 
A  +  B2  S[(72  (x  —  D)],  x>D. 


(25) 


Although  this  form  is  written  with  six  free  parameters,  in  actuality,  to  ensure  continuity  of  the  first 
derivate  at  the  transition  point,  one  degree  of  freedom  is  absorbed,  whereby  B\  C\  =  B2C2. 


To  model  the  baseline  value  we  use, 


V0{Q)  =  Ss[log2(Q),  0.870,  0.370,  0.085,  0.355,  1.545,  -1.00].  (26) 

Using  this  equation,  the  computed  V  results  based  on  the  integration  technique  described  in  the 
previous  section  were  rendered  into  Va  by  subtracting  Vo(Q).  Next,  it  was  recognized  that  the 
different  families  of  Va  results  could  be  better  compared  by  normalizing  each  set  of  results  at 
different  Q  values  by  dividing  through  by  the  maximum  value  of  Va  in  each  subset.  This 
maximum  value  invariably  occurred  at  the  maximum  X  value  used  (X  =  103)  at  the  peak  point 
u>p(Q,  X).  This  maximum  Va  (Vm)  is  given  approximately  by  the  relation, 

Vm(Q)  =  S[log2(Q),  0.213,  0.072,  1.525,  0.150],  (27) 

Both  Vq(Q)  and  Vm{Q )  are  plotted  in  figure  3  versus  database  derived  results. 
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iog2(Q) 


Figure  3.  Plots  of  Vq(Q)  and  Vm{Q)- 


Division  by  Vm(Q)  produces  sets  of  similar  functions  at  each  Q,  Va  =  Va/Vm-  Figures  4  and  5 
compare  two  sets  of  such  functions. 


UJ 

Figure  4.  Plots  of  Va {Q,  X,  oS)  for  Q=l/16. 
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Figure  5.  Plots  of  Va (Q,  A,  u)  for  Q=16. 


Both  plots  illustrate  sets  of  curves  of  similar  shape.  Both  figures  include  a  plot  of  the  path 
corresponding  to  the  position  and  height  of  the  peak  of  each  curve  as  a  function  of  increasing  X 
for  that  value  of  Q.  Apparently  the  tracks  of  these  peak  curves  are  nearly  identical,  but  different 
points  along  the  paths  are  associated  with  different  X  values.  As  an  illustration,  figure  6  plots  the 
u  components  of  these  curves  for  various  Q  values  as  functions  of  x  =  log10(A")  [uP(Q,  A")]. 

To  model  the  behavior  of  uP,  the  behavior  of  uP(Q  =  1,  A)  may  be  used  as  the  main  guide 
function  and  shifted  in  x  =  log10(A)  by  a  function  dependent  on  0  alone: 


u)P(Q,X)  =  uP [Q  =  1,  X/XU(Q)\.  (28) 

uP(Q  =  1,  A)  =  Ss[log10(A),  0.360,  -0.133,  -0.322,  3.450,  1.424,  1.320],  (29) 

log10[Aw(Q)]  =  S[log2(Q),  -0.004,  0.091,  1.750,  -0.05],  (30) 

The  resulting  curve  fit  for  uP(Q ,  A")  is  shown  for  the  Q  —  1  case  in  figure  7.  Several  shifted 
versions  of  other  data  sets  based  on  the  A S(Q)  curve  are  also  plotted. 
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logio(X) 


Figure  6.  Plots  of  lop(Q,  X )  for  various  Q  values. 


logio(X) 


Figure  7.  Plots  of  c op(Q  =  1,  A")  for  various  x-shifted  curves  based  on  the 
up{Q  =  1,  X)  parametric  curve. 
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Next,  consider  the  vertical  component  of  the  peak  curve.  Define  this  function  as 
Vp  (■ Q ,  A")  =  Va  (Q,  X ,  ujp) .  The  behavior  of  Vp  {().  X )  is  illustrated  for  several  values  of  Q  in 
figure  8.  The  salient  features  of  these  plots  are  that  there  is  a  similar  shift  in  x  =  log10(A) 
between  the  different  curves  along  with  a  vertical  scale  factor. 


logio(X) 


Figure  8.  Plots  of  Vp(Q.  X )  for  various  values  of  Q. 


To  model  this  behavior  let  us  again  characterize  Vp(Q,  X)  by  first  considering  the  peak  function 
for  Q  —  1: 


Vx{x)  =  VP(Q  =  l,x)  =  Ss(x,  0.122,  0.044,  0.082,  4.257,  2.300,  1.167).  (31) 

Using  this  baseline  function  we  rescale  and  shift  it  according  to: 

Vp(q,  x)  =  AR(q)  Vi  [x  -  Xp(q )],  q  =  log 2(Q),  x  =  log10(X);  (32) 

AR(q)  =  S(q,  +1.051,  +0.3565,  +1.600,  +0.150);  (33) 

xP(q)  =  S(q,  +0.005,  -0.0880,  +1.550,  -0.150).  (34) 
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Given  the  peak  value  function,  Vp,  and  the  peak  location,  uP,  one  may  then  approximately  model 
the  functional  dependence  of  V(Q,  X,  tu).  However,  the  shape  of  the  resulting  function  must  still 
be  parameterized  by  indexing  this  shape  relative  to  the  normalized  peak  value, 

UP(Q,X)  =  Vp(Q,X)/VP(Q,  1000). 

To  model  the  shape  of  these  functions,  a  scaled  set  of  Legendre  polynomial  functions  was 
employed.  Let  Po(u)  =  1,  P\{u)  =  u,  and  use  equation  5.5.1  of  Press  et  al.  (5): 

(n  +  1)  Pn+i(u)  =  (2 n  +  1  )uPn(u)  -  nPn-i(u).  (35) 

The  u  variable  of  Legendre  polynomials  ranges  between  -1  and  +1,  where  it  is  evenly  weighted. 
However,  because  0  <  uj  <  1,  we  need  to  introduce  new  functions, 

Ln{oj)  =  v/2^TlP„(2w  -  1),  (36) 


such  that, 

/  Ln(uJ )  dcU  $mni  (37) 

Jo 

expresses  the  orthonormal  nature  of  the  resulting  rescaled  functions. 

In  addition,  the  function  expansion  in  weighted  Lrn  functions  may  be  reduced  to  fewer  terms  if 
the  interval  (0...1)  in  u  is  divided  into  two  subsections  (CL.cup)  and  (cup...l).  Then,  the 
expansion  of  Va(Q,  X,  uj)  may  be  expressed  as, 

n=N  r  'j 

Va(Q,X,u>)  «  n[Up(Q,X)]L2n\-  — ,  u<up- 

n= 0  [ZUJp^QjX  JJ 

Va{Q,X,u)  «  nY^D2n[UP(Q,X)]L2n{l-  [1~^  V  u  >  up.  (38) 

n=0  l  2  [1  -  Up[Q,  X)\  ) 

This  division  is  significant  because  only  even  order  scaled-Legendre  functions  are  needed,  and 
also  at  u  >  ujp  fewer  significant  terms  are  needed.  In  particular,  the  series  has  been  truncated  at 
N  =  6.  The  constants  GV,,  are  then  produced  by  integrating  the  scaled  functions  multiplied  by 
each  of  the  rescaled  Legendre  functions. 

Though  some  scatter  appears  when  performing  these  calculations,  the  similarities  between  the 
functions  at  different  Q  values  are  evident.  These  results  are  plotted  in  figures  9  and  10  for  the 
C2n  and  D-2n  coefficients. 
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Up 


Figure  9.  Plots  of  Cm  =  Cin  expansion  coefficients  and  parametric  curve 
fits  to  the  results. 


Up 


Figure  10.  Plots  of  Dm  =  Dm  expansion  coefficients  and  parametric 
curve  fits  to  the  results. 
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For  the  C^n  plots,  several  coefficients  appear  to  exhibit  hyperbolic  dependence,  while  of  the  D2n 
coefficients  only  D0  appears  to  exhibit  hyperbolic  behavior.  These  hyperbolic  forms  were 
modeled  using  the  fitting  function, 

Cm  =  &m  (Up  Cm)  T  bm  \J (Up  £m T  cm  T  dm.  (39) 


The  curves  plotted  in  figures  9  and  10  are  based  on  this  equation.  Coefficients  used  in  these 
equations  are  detailed  in  tables  1  and  2.  Table  1  details  the  coefficients  associated  with  the 
hyperbolic  cases.  The  remaining  coefficients  are  considered  linear  in  Up,  where  the  same  model 
terms  are  used,  but  where  constants  bm,  cm,  and  em  are  all  set  to  zero. 


Table  1.  Hyperbolic  Model  Cm  and  Dm  Coefficients. 


Coefficient 

dm 

bm 

Cm 

dm 

Cm 

Co 

20.1020 

19.5598 

0.0028 

0.5700 

1.0500 

c2 

2.1140 

2.4810 

0.0041 

-0.3503 

0.9997 

c4 

-0.3391 

-0.3361 

0.0050 

0.0051 

0.8141 

c6 

-0.3533 

-0.3721 

0.0000 

0.0131 

0.9201 

c8 

-0.3675 

-0.3987 

0.0001 

0.0155 

0.9632 

D0 

6.4460 

5.9010 

0.0006 

0.5810 

0.9999 

Table  2.  Linear  Model  Cm  and  Dm 
Coefficients. 


Coefficient 

dm 

dm 

Cio 

0.01013 

-0.00448 

Cl2 

0.00558 

-0.00256 

d2 

-0.30465 

-0.00055 

Da 

0.03405 

-0.00396 

De 

-0.01269 

0.00486 

d8 

0.00429 

-0.00153 

Dio 

-0.00444 

0.00191 

D\2 

0.00055 

-0.00016 

The  root  mean  square  (RMS)  error  resulting  from  the  use  of  these  coefficients  for  V  results 
computed  at  9  Q  values,  41  X  values,  and  105  u  values  was  0.006795.  This  result  was  further 
improved  by  passing  the  resulting  model  into  an  algorithm  utilizing  a  Stochastic  Parallel  Gradient 
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Descent  (SPGD)  approach  of  Vorontsov  and  Sivokon  (6).  The  result  of  running  this  algorithm 
tended  to  optimize  the  coefficients,  reducing  the  overall  RMS  error  described  above  to  0.002728. 
Based  on  these  results  we  must  redefine  the  various  equations  previously  presented.  These  are 


expressed  using, 

V0(Q)  =  Ss[log2(Q),  0.8785,  0.3726,  0.0849,  0.3733,  1.5700,  -0.9844];  (40) 

Vm{Q)  =  S[log2(Q),  0.2205,  0.0696,  1.5387,  0.1502];  (41) 

uP(Q  =  l,X)  =  Ss[log10(X),  0.3506,  -0.1332,  -0.3177,  3.4505,  1.4065,  1.2738];  (42) 

log10[Xw(Q)]  =  5[log2(Q).,  0.0040,  -0.0877,  1.7450,  -0.0502];  (43) 

Vx{x)  =  VP(Q  =  l,x)  =  Ss{x,  0.1216,  0.0438,  0.0823,  4.2957,  2.3225,  1.1787);  (44) 

AR(q)  =  S(q,  +1.0657,  +0.3521,  +1.5964,  +0.1498);  (45) 

xP(q)  =  S(q,  -0.0050,  +0.0893,  +1.5499,  -0.1500).  (46) 


With  these  changes,  the  results  in  tables  1  and  2  are  replaced  by  the  coefficients  of  tables  3  and  4. 
The  primary  difference  between  these  coefficients  and  the  original  parameter  values  appears  to  be 
due  to  the  restrictions  on  the  range  of  the  Q  and  X  parameters  to  the  intervals  1/16  -  16  and  OH¬ 
IO3,  respectively.  The  SPGD  method  appears  quite  capable  of  optimizing  the  fitting  properties  if 
presented  with  a  coefficient  set  that  is  fairly  close  to  the  optimal  settings.  On  the  other  hand,  if 
the  SPGD  method  is  supplied  with  an  initial  coefficient  set  that  is  initially  a  rather  poor  fit  to  the 
available  data,  then  SPGD  does  very  poorly  in  improving  the  fit  to  the  data. 


Table  3.  SPGD  Modified  Hyperbolic  Model  Cm  and  Dm 
Coefficients. 


Coefficient 

&m 

bm 

Cm 

dm 

Cm 

Go 

20.2124 

19.6324 

0.0031 

0.5751 

1.0686 

g2 

2.1118 

2.4928 

0.0042 

-0.3515 

0.9808 

g4 

-0.3413 

-0.3355 

0.0051 

0.0051 

0.7996 

g6 

-0.3529 

-0.3769 

0.0000 

0.0131 

0.9041 

g8 

-0.3696 

-0.3982 

0.0000 

0.0153 

0.9504 

Do 

6.4372 

5.8919 

0.0005 

0.5788 

1.0341 
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Table  4.  SPGD  Modified  Linear 
Model  Cm  and  Dm 
Coefficients. 


Coefficient 

dm 

dm 

C'io 

0.0 1 009 

-0.00452 

Cl  2 

0.00550 

-0.00257 

D2 

-0.30273 

-0.00057 

da 

0.03335 

-0.00413 

d6 

-0.01271 

0.00485 

d8 

0.00425 

-0.00154 

Dio 

-0.00452 

0.00186 

D\2 

0.00055 

-0.00019 

4.  Results 


The  new  model  of  V(Q,  X,  u)  allows  us  to  compare  the  new  results  with  the  previous  model 
results.  Perhaps  the  most  telling  examples  are  comparisons  of  the  same  cases  presented  in 
figure  6  of  Tofsted  ( 1 ),  reproduced  here  as  figure  11.  The  new  results  are  plotted  in  figure  12 
below.  In  the  prequel,  several  cases  showed  “super-resolution”  characteristics  by  exceeding  unity 
at  u  values  greater  than  0.7.  The  new  curves  never  violate  the  diffraction  limit,  as  argued  by 
Charnotskii  (4),  but  also  produce  increased  values  in  the  functions  at  intermediate  frequencies 
0.3  <  u  <  0.7  because  V (Q,  X,  u)  is  larger  at  intermediate  lu  values  than  it  is  at  the  knee  of  the 
curves. 

The  difference  between  the  prequel  paper’s  super-resolution  results  and  the  diffraction  limited 
results  obtained  here  are  also  illustrated  in  figure  13  for  a  presumed  worst  case  fit  of  the  previous 
analysis.  That  scenario  is  Q  =  4  and  X  =  2.  This  combination  produced  Msa{ w)  >  1.25  at 
lu  >  0.95  in  the  previous  analysis.  As  the  figure  shows,  the  prequel’s  V  =  1.0395  produces  a 
function  closer  to  the  true  MTF  than  the  a  =  1.0  Fried  case  up  until  uj  =  0.8.  The  new  function 
is  virtually  identical  to  the  calculated  curve.  Hence,  the  prequel’s  functional  form  still  provides  a 
useful  approximation,  even  though  it  violates  the  diffraction  limit  at  extremely  high  frequency 
(u  >  0.87). 
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(a)  (b) 


(c)  (d) 


Figure  11.  Atmospheric  MTF’s  plotted  for  Q  cases  (a)  1/2,  (b)  1,  (c)  2,  and 
(d)  4,  respectively,  from  the  prequel  paper  where 
super-resolution  conditions  were  predicted  where  results 
exceeded  1.0. 


The  various  models  can  be  compared  against  the  complete  set  of  short-exposure  computed  MTF 
values,  i.e.,  the  calculated  V ( Q ,  X ,  ui)  data  set.  To  produce  a  weighted  sum  squared  error,  the 
sum  of  the  weighted  squared  difference  between  an  estimate  and  the  computed  result,  u  5M’2,  can 
be  compared  for  five  different  methods  of  estimating  the  MTF.  These  estimation  techniques  are: 
(1)  Fried’s  far- field  case  (a  =  1/2),  (2)  Fried’s  near- field  case  ( a  =  1),  (3)  a  Fried  best-fit  case 
(based  on  either  a  =  1/2  or  1,  depending  on  which  case  best  fits  the  data),  (4)  the  prequel’s 
analytical  model  V ( Q ,  X),  and  (5)  the  current  paper’s  extended  analytical  V (Q.  X,  u)  model. 

For  each  case  the  error  between  the  computed  MTF  and  the  estimated  result  is  computed  (call  it 
5M).  The  root  mean  squared  error  was  then  normalized  according  to, 

Rc  =  1^2  5M?  Ui  A Ui/  W  Acoi,  (47) 

V  i  i 

where  C  is  the  case  number.  The  resulting  five  cases  produced  RMS  errors  Ar,  =  0.01409, 

N2  =  0.03518,  IV3  =  0.01079,  IV4  =  0.00643,  and  N5  =  0.00022.  Clearly  the  worst  case  is  the 
Fried  a  =  1/2  far- field  case.  Of  course,  the  cases  studied  focused  on  the  transition  range  from 
small  to  large  Q  values.  If  Q  were  extremely  small,  the  a  =  1/2  case  might  be  the  best  choice. 
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Figure  12.  Atmospheric  MTF’s  plotted  for  Q  cases  (a)  1/2,  (b)  1,  (c)  2,  and 
(d)  4,  respectively.  Super-resolution  cases  of  the  prequel  model 
are  eliminated  due  to  high  frequency  decay. 


However,  even  when  permitting  a  choice  between  the  near-held  and  far-held  cases,  the  Fried 
model  exhibits  signihcantly  greater  error  than  the  case  4  results  using  the  prequel’s  V {().  X ) 
function.  And  of  course,  the  case  5  model  is  obviously  virtually  equivalent  to  the  full 
multi-dimensional  integral  of  equations  13  and  14,  being  ~  30  times  better  than  the  prequel’s 
performance. 


5.  Conclusions 


The  analysis  performed  in  this  paper  completes  a  research  arc  beginning  nearly  a  decade  ago  with 
the  publication  of  Tofsted  (7).  For  many  years  researchers  (c.f.,  MacDonald  and  Cain  (VS’),  and 
Dunphy  and  Kerr  (9))  have  alluded  to  weaknesses  in  Fried’s  (2)  short-exposure  model.  The  path 
of  this  arc  led  through  Tofsted  (/)  which  highlighted  the  missing  term  in  Fried’s  approach.  This 
paper  completes  the  arc,  following  Charnotskii’s  (4)  observation  that  the  theory  cannot  support 
super-resolution,  and  resulted  in  a  new  understanding  of  the  angular  frequency  in  the  V  function. 
As  before,  an  analytical  expression  was  derived  to  describe  the  new  dependence,  though  involving 
considerably  more  detail.  In  the  process  the  degree  of  accuracy  of  the  new  method  was 
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Figure  13.  Atmospheric  MTF’s  plotted  for  Q=4  and  X=2  comparing 

Fried’s  model  (red  dashed)  for  <alpha>=l,  Tofsted  approximate 
model  (blue  dashed)  for  V=  1.0395,  and  full  model  (black  solid) 
results,  along  with  system  MTF  (dashed  green)  curve.. 


compared,  along  with  the  prequel  version  and  Fried’s  original  analysis  with  the  more  extensive 
calculations,  revealing  that  the  new  method  has  virtually  no  error,  while  the  prequel  method  has 
one-third  the  error  of  Fried’s  original  method  using  a  —  1. 

These  new  results  can  also  be  directly  coupled  to  the  results  of  a  recent  paper  (10)  that  studied  the 
influence  of  path- varying  turbulence  strength  for  two  classes  of  path  geometry:  slant  path  over 
flat  terrain  and  flat  path  over  valley  terrain.  In  both  cases  it  was  found  that  changes  in  the 
behavior  of  the  wave  and  phase  structure  functions  with  path  varying  turbulence  could  be 
parameterized  through  scaling  functions  based  on  path  geometry  to  modify  the  X  and  Q 
parameters  of  these  structure  functions.  Since  the  results  of  this  paper  continue  to  use  the  same 
structure  functions,  the  same  techniques  may  be  applied  to  adjust  the  X  and  Q  parameters  of  the 
new  MTF  model  to  account  for  path- varying  turbulence. 

The  key  missing  element  of  the  current  and  previous  studies  has  been  and  is  the  lack  of 
consideration  of  inner  scale  influences  on  the  structure  functions.  However,  since  the  focus  of 
these  studies  has  been  the  impacts  of  turbulence  on  aided  optical  systems,  in  most  cases  inner 
scale  effects  are  minor.  In  contrast,  outer  scale  effects  may  be  significant  on  the  angle-of-arrival 
variance,  but  only  for  larger  aperture  optics.  For  systems  considered  here  that  are  less  than 
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10-inches  in  aperture  diameter,  the  effects  of  outer  scale  on  angle-of-arrival  variance  will  be  less 
than  10%  (see  e.g.,  Tofsted  (10),  which  is  typically  less  (under  most  conditions)  than  the 
uncertainty  in  the  measurement  of  C*  itself. 
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